Pipeline para el análisis de datos de transcriptómica espacial GeoMx DSP
Carme Delicado Mercader
2025-12-04
El desarrollo de este pipeline está basado en las dos herramientas existentes para analizar datos provinentes de la plataforma GeoMx DSP: GeoMxTools y StandR.
Se hace uso de GeoMxTools, el pipeline desarrollado por Nanostring que trabaja con NanoStringGeoMxSet como objeto principal, para la carga y visualización de los datos.
Posteriormente se implementa el resto del pipeline con StandR, que fue desarrollado en 2023, y contiene métodos más robustos. Trabaja con SpatialExperiment como objeto, lo que proporciona facilidad a la hora de integrar funciones de otros pipelines para los análisis posteriores.
Las viñetas de ambos pipelines se detallan aquí:
1. Preparación del entorno de trabajo
1.1 Instalación de paquetes necesarios
Para la ejecución de este Pipeline se hará uso de funciones del pipeline de GeomxTools y de StandR.
Cargamos las librerías necesarias.
1.2 Cargar librerías
#Librerías específicas de StandR
library(standR)
library(SpatialExperiment)
#Librerías específicas de GeomxTools
library(GeomxTools)
#Librerías para ordenar datos y encadenar funciones
library(tidyverse)
library(dplyr)
#Librerías para graficar
library(ggplot2)
library(networkD3)
library(ggforce)
library(ggalluvial)
library(ggrepel)
#Descomprime archivos gzip (.gz) y .tar
library(R.utils)
#Cálculo logCPM
library(edgeR)
#tablas
library(knitr)
#Tabla interactiva
library(DT)
#Excel
library(readxl)
library(writexl)
#Limma pipeline
library(limma)2. Dataset
2.1 Introducción al dataset
El GBM es el tumor cerebral primario maligno más agresivo y común entre la población adulta. Tiene un crecimiento muy rápido y presenta resistencia a las terapias estándares, por lo que en la mayoría de casos los pacientes experimentan recurrencia del tumor.
La inmunoterapia ha supuesto un avance significativo en el tratamiento de varios tipos de cáncer, impulsando respuestas al aprovechar y potenciar los mecanismos naturales del sistema inmunitario para reconocer y eliminar las células tumorales. Un ejemplo de inmunoterapia para tratar los casos recurrentes de GBM es el bloqueo entre el PD-1 y su ligando PD-1L. Cuando se bloquea esta interacción, los inhibidores de PD-1 restauran la actividad de las células T, permitiendo que las células T citotóxicas eliminen las células tumorales.
El éxito de esta práctica es confuso en cuanto a resultados, hasta la fecha no se ha demostrado si el bloqueo de PD-1 induce cambios significativos en el microambiente tumoral (TME) del GBM, y por lo tanto, tiene un impacto en la supervivencia de los pacientes con GBM.
Para abordar esta cuestión se han llevado a cabo varios estudios recientemente, en los cuales los pacientes reciben tratamiento con un neoadyuvante (anti-PD-1) antes de la resección tumoral, lo que permite analizar el tejido tumoral poco después del tratamiento. De esta forma, se intenta evaluar si ha habido una activación por parte del sistema inmunitario tras el bloqueo del ligando.
Por lo tanto, el objetivo del estudio es utilizar la técnica GeoMx DSP para evaluar el efecto del bloqueo de PD-1 en células tumorales y Tumor Associated Macrophages (TAM). El análisis de estas zonas, corresponde con la aplicación más común de GeoMx DSP, el estudio de la población tumoral y su microambiente.
En este estudio, se seleccionan 30 pacientes con GBM recurrente.
20 pacientes reciben el tratamiento con Nivolumab, el neoadyuvante, que bloquea la interacción de PD-1 y su ligando.
10 pacientes son controles emparejados no tratados.
Los 20 pacientes, antes de ser seleccionados, han sido testeados para comprobar la participación y activación de células T en cuanto se les bloqueaba la interacción PD-1 - PD-1L.
Se pretende observar si la actividad inmunitaria de las células T va acompañada de cambios transcripcionales en las células tumorales y los TAM. Se seleccionan zonas con alta densidad de SOX2+, que son indicadores de células tumorales, y IBA1+, indicadores de TAM.
Los hallazgos publicados por el estudio, de donde se extrae el dataset, dictan que no se observan cambios transcripcionales significativos en las células tumorales o los TAM tras el bloqueo de PD-1. Por lo tanto, sugieren que el impacto de la monoterapia neoadyuvante con PD-1 en el GBM recurrente puede ser limitado. Para más información: https://pubmed.ncbi.nlm.nih.gov/40960500/
2.2 Cargar dataset
Se accede a los datos mediante el siguiente URL: https://zenodo.org/records/16839828
Se encuentran diferentes archivos disponibles para descargar:
RAW.tar.gz: archivos dcc (GeoMx output) de todas las AOI (i.e., raw data)
metadata.csv.tar.gz: archivo separado por comas de las AOIs utilizados en ester análisis.
sampleAnno.xlsx.tar.gz: archivo excel que contiene la anotación del fichero dcc.
Hs_R_NGS_WTA_v1.0.pkc.tar.gz: archivo pkc que contiene la asociación de sonda-gen.
normalized.txt.tar.gz: archivo separado por tabulaciones que contiene la matriz de conteos normalizada y corregida por lotes de las AOIs que se usaron en este análisis.
El objetivo de este experimento es implementar el pipeline completo, que cuenta con el análisis de datos desde el preprocesamiento hasta los análisis posteriores. Por dicha razón, es de nuestro interés descargarse los datos curdos y no los normalizados, a parte de todo el resto de archivos.
Se configura el directorio de trabajo.
## [1] "Hs_R_NGS_WTA_v1.0.pkc.tar.gz" "metadata.csv.tar.gz"
## [3] "RAW.tar.gz" "sampleAnno.xlsx.tar.gz"Los archivos tienen extensiones .tar y .gz.
.tar sirve para empaquetar varios archivos en uno solo.
.gz comprime el archivo único usando gzip.
Para acceder a los datos se deben descomprimir los archivos mediante la función untar().
untar(file.path(datadir,"RAW.tar.gz"), exdir = file.path(datadir,"RAW_data"))
untar(file.path(datadir,"metadata.csv.tar.gz"), exdir = file.path(datadir,"metadata"))
untar(file.path(datadir,"sampleAnno.xlsx.tar.gz"), exdir = file.path(datadir,"sampleAnno"))
untar(file.path(datadir,"Hs_R_NGS_WTA_v1.0.pkc.tar.gz"), exdir= file.path(datadir,"PKCFiles"))Se listan los archivos de las carpetas descomprimidas para conocer su contenido.
## [1] "DCC-20230619"## [1] "metadata.csv"## [1] "sampleAnno.xlsx"## [1] "Hs_R_NGS_WTA_v1.0.pkc"Los archivos de trabajo se obtienen directamente en tres de los
cuatro casos. Sin embargo, en la carpeta denominada
DCC-20230619, se encuentran los archivos con extensiones
.dcc correspondientes a cada AOI.
Se listan los cinco primeros archivos de DCC-20230619,
para conocer su interior.
## [1] "DSP-1001660012264-D-A01.dcc" "DSP-1001660012264-D-A02.dcc"
## [3] "DSP-1001660012264-D-A03.dcc" "DSP-1001660012264-D-A04.dcc"
## [5] "DSP-1001660012264-D-A05.dcc"Para cargar los datos el pipeline de GeoMxTools utiliza la función
readNanoStringGeoMxSet() que unifica los archivos en un
único objeto y que necesita los siguientes argumentos:
dccFiles <- Un vector de caracteres que contiene la
ruta al fichero DCC, que contiene los archivos .dcc.
pkcFiles <- Un vector de caracteres que representa la
ruta correspondiente al fichero PKC .pkc.
phenoData <- Un dataframe que contiene los datos
fenotípicos. Si éste y phenoDataFile están disponibles, se
usa phenoData.
phenoDataFile <- Una cadena de texto que representa
la ruta al archivo correspondiente de datos fenotípicos. El archivo de
datos puede ser una hoja de laboratorio, un archivo Excel u otro archivo
delimitado. Se recomienda usar la hoja de laboratorio en el mismo orden
exacto en el que se proporcionan las muestras.
phenoDataDccColName <- Cadena de texto que identifica
la columna del archivo phenoData o
phenoDataFile que contiene el identificador único de cada
muestra.
Se crean estos objetos.
dccFiles <- dir(file.path(datadir, "RAW_data/DCC-20230619"),
pattern = ".dcc$", #guarda solo docs acabados con .dcc
full.names = TRUE,
recursive = TRUE)
pkcFiles <- file.path(datadir,"PKCFiles", "Hs_R_NGS_WTA_v1.0.pkc")
phenoDataFile <- file.path(datadir,"sampleAnno","phenoDataFile.xlsx")El fichero Excel viene con 21 filas que se deben eliminar antes de
crear el objeto datos.
sampleanno_excel <- read_excel(path = paste(datadir,"sampleAnno","sampleAnno.xlsx", sep = "/"),
skip = 21 )
write_xlsx(sampleanno_excel, path = paste(datadir,"sampleAnno","phenoDataFile.xlsx", sep = "/"))Se unifican los objetos en un único objeto llamado datos
de clase NanoStringGeoMxSet.
datos <- readNanoStringGeoMxSet(dccFiles = dccFiles,
pkcFiles = pkcFiles,
phenoDataFile = phenoDataFile,
phenoDataSheet = "Sheet1",
phenoDataDccColName = "Sample_ID")## Warning in readNanoStringGeoMxSet(dccFiles = dccFiles, pkcFiles = pkcFiles, :
## The following DCC files had no counts: DSP-1001660017941-E-H03.dcc,
## DSP-1001660017943-G-G04.dcc, DSP-1001660017943-G-H02.dcc, These will be
## excluded from the GeoMxSet object.## Warning in readNanoStringGeoMxSet(dccFiles = dccFiles, pkcFiles = pkcFiles, :
## DCC files missing for the following: TOTAL AREA PLATE A.dcc, TOTAL AREA PLATE
## B.dcc, TOTAL AREA PLATE C.dcc, TOTAL AREA PLATE D.dcc, TOTAL AREA PLATE E.dcc,
## TOTAL AREA PLATE F.dcc, TOTAL AREA PLATE G.dcc, TOTAL AREA PLATE H.dcc, These
## will be excluded from the GeoMxSet object.Aparecen avisos que indican la eliminación de aquellos ficheros que no contienen la información necesaria para ser procesados e incluidos en el estudio.
2.3 Ajustes y visualización del dataset
Se accede a la información básica del conjunto de datos.
El objeto datos es un objeto de clase
NanoStringGeoMxSet, que proviene del paquete
GeomxTools.
Para poder explorarlo hay que conocer las funciones específicas del
pipeline GeomxTools.
Mediante la función exprs() se puede acceder a la matriz
de conteos, que representa el número total de sondas que contiene cada
AOI.
## [1] 18815 713Se cuenta con un dataset que contiene 18815 sondas y 713 AOIs.
Para acceder a las anotaciones de los genes se utiliza la función
fData(datos), que establece la relación sonda-gen.
Las Áreas De Iluminación (AOI) se conocen también como segmentos. Una vez elegidas las Regiones De Interés (ROI) del tejido, estas se pueden segmentar en las distintas áreas celulares que se quieren secuenciar. Solo estas áreas se van a iluminar con UV para que se produzca la liberación del identificador único del mRNA, lo que nos aportará la información sobre los mRNA se encuentran en esa área. De esta forma, una ROI puede contener más de una AOI, o bien, si una ROI solo tiene una AOI, los conteos de esa ROI coinciden con los de la AOI.
Para acceder a las anotaciones de las AOIs se usa la función
pData().
Se exploran las distintas AOI de cada tipo, es decir, las zonas celulares que han sido iluminadas.
Además, se analiza el número de ROIs totales y la cantidad de AOIs presentes en cada ROI.
##
## GFAP-aoi-001 Iba1-aoi-001 Sox2-aoi-001
## 224 258 223## [1] 205##
## C01R1is C01R2g C01R2is C01R3g C01R3ig C01R4g C01R4ig C02R1is C02R2g C03R1is
## 10 3 4 1 2 1 2 6 3 10
## C03R2g C03R2is C03R3g C03R4g C04R1is C04R2g C04R2is C04R3ig C04R4ig C04R5g
## 4 2 1 1 10 4 2 2 2 1
## C05R1g C05R1is C05R2g C05R2is C05R3is C05R4g C05R5g C05R6g C06R1is C06R2g
## 1 10 4 4 4 1 1 1 8 3
## C06R2is C06R3g C06R4g C06R5g C07R1is C07R2g C07R2ig C07R2is C07R3g C07R3is
## 2 1 1 1 10 3 2 2 1 2
## C07R4g C08R1is C08R2g C08R2is C08R3g C08R3is C08R4g C09R1is C09R2g C09R2ig
## 2 8 3 2 1 2 2 8 3 2
## C09R3ig C10R1is C10R2g C10R2ig C10R2is C10R3is C10R4is C10R5g C10R6g C10R7g
## 2 10 3 2 2 2 2 1 1 1
## C10R8g C11R1is C11R2g C11R2ig C11R2is C11R3g C11R4g C12R1is C12R2g C12R2ig
## 1 10 2 2 4 2 1 8 2 2
## C12R2is C12R3g C12R4g C13R1is C13R2g C13R2is C13R3g C13R3is C13R4g C14R1g
## 2 1 1 8 3 2 1 2 1 1
## C14R1is C14R2g C14R2is C14R3g C14R4g C15R1is C15R2g C15R2is C15R3g C15R4g
## 6 2 2 1 1 10 4 2 1 1
## C16R1is C16R2g C17R1is C17R2g C18R1is C18R2g C19R1is C19R2g C20R1is C20R2g
## 8 4 6 3 8 4 10 5 8 3
## C20R2ig C21R1is C21R2g C21R2ig C21R2is C21R3g C21R4g C22R1is C22R2g C22R2is
## 2 10 2 4 2 1 1 10 4 2
## C22R3g C23R1is C23R2g C23R2ig C23R3is C23R4g C24R1is C24R2g C24R3g C25R1is
## 1 10 3 2 2 1 8 4 1 12
## C25R2g C25R2is C25R3g C25R4g C26R1is C26R2g C26R2is C26R3g C26R3ig C26R4ig
## 3 3 2 2 10 4 2 1 2 2
## C27R1is C27R2g C27R2is C27R3g C28R1g C28R1is C28R2g C28R2ig C28R3is C28R4g
## 10 4 2 1 1 8 3 2 2 1
## C28R5ig C29R1is C29R2g C29R2ig C29R2is C29R3g C29R4g C30R1is C30R2g C30R2ig
## 2 8 1 4 2 1 1 12 2 4
## C30R2is C30R3g C30R4g C31R1is C31R2g C31R2is C31R3g C31R4g C32R1is C32R2g
## 4 2 2 10 4 2 1 1 8 3
## C32R2ig C33R1is C33R2g C33R2ig C33R2is C33R3g C33R4g C34R1g C34R1is C34R2g
## 1 12 1 4 4 2 1 1 8 2
## C34R2ig C34R3g C35R1is C35R2g C35R2ig C35R3g C36R1is C36R2g C36R2ig C37R1is
## 4 1 10 2 6 1 10 2 4 8
## C37R2g C37R2ig C38R1is C38R2g C38R2ig C38R2is C38R3is C38R4is C38R5g C38R6g
## 2 4 10 3 2 2 2 2 1 1
## C38R7g C38R8g C39R1g C39R1is C39R2g C39R2ig C39R2is C39R3g C39R4g C40R1g
## 1 1 1 8 3 2 2 1 1 1
## C40R1is C40R2g C40R2is C40R3is
## 6 3 2 2Las AOI se dividen en 224 GFAP, 258 Iba1 y 223 Sox2, que provienen de un total de 205 ROIs.
En el apartado 2.2 se ha descargado el archivo metadata,
no utilizado hasta ahora. Este contiene las anotaciones de las AOIs
utilizadas para el análisis de este experimento. Se descarga y se
explora para conocer qué diferencias frente al
annotation_AOI contiene.
## [1] 482 63metadata contiene 482 AOI a diferencia de las 713
anteriores.
También contiene 63 variables a diferencia de las 17 anteriores.
Se exploran las distintas AOI de cada tipo y el número de ROIs totales.
##
## Iba1 Sox2
## 259 223## [1] 103Las AOI de metadata se dividen en 259 Iba1 y 223
Sox2.
Se concluye que el experimento no agrega las 224 de tipo GFAP que se
encuentran en annotation_AOI. Se reducen también el número
total de ROIs a 103.
Se obtiene un AOI más de tipo Iba1 en comparación con el resultado
obtenido en annotation_AOI.
Se crea una tabla de los metadata para tener un
conocimiento de la información que recoge.
metadata <- metadata %>% mutate(across(where(is.character),as.factor))
library(DT)
datatable(
metadata,
extensions = "Scroller",
options = list(
deferRender = TRUE,
scrollX = TRUE,
scrollY = 500,
scroller = TRUE,
pageLength = 10
),
class = "compact stripe"
)Para empezar con el análisis se añaden las variables faltantes de
metadata al objeto datos, al igual que se
eliminan las AOI de GFAP del objeto datos.
datos_con_GFAP <- datos
anno_AOI_filtrado <- annotation_AOI %>%
filter(Segment %in% c("Sox2", "Iba1"))
selected_cols <- rownames(anno_AOI_filtrado)
datos <- datos_con_GFAP[, selected_cols]
dim(datos)## Features Samples
## 18815 481Se obtienen 481 AOI en datos.
Se crea una columna que será la que se comparará entre los objetos
datos y metadata.
(Recordatorio: el objeto annotation_AOI sale de
pData(datos)).
pData(datos)$unique_ID <- paste(
pData(datos)$`Slide Name`,
pData(datos)$Roi,
pData(datos)$Segment,
sep = "_"
)
anno_AOI_filtrado <- pData(datos)Se enfrenta metadata$unique_ID y
pData(datos)$unique_ID para obtener la AOI sobrante.
## [1] "SaraWTATMA4_C25R2is_Iba1"Se elimina esta AOI de metadata.
## [1] 481 63Se juntan pData(datos) y metadata.
Es importante quedarse con la nomenclatura de las AOI de las filas de
pData(datos), ya que se corresponden con el mismo nombre
que tienen las columnas de los conteos exprs(datos).
anno_AOI_filtrado <- anno_AOI_filtrado %>%
rownames_to_column(var = "AOI_name")
Annotation_AOI_completo <- anno_AOI_filtrado %>%
left_join(metadata, by = "unique_ID")
rownames(Annotation_AOI_completo) <- Annotation_AOI_completo$AOI_name
Annotation_AOI_completo$AOI_name <- NULL
#verificación
all.equal(rownames(Annotation_AOI_completo), colnames(exprs(datos)))## [1] TRUE## [1] 481 80Se grafican las réplicas de cada anotación para comprobar que se obtiene un mínimo de entre 3-5 réplicas para cada grupo.
Grupo = SegmentLabel + trial_setting + tumor_setting
pData(datos)$Grupo_gmx <- paste(pData(datos)$SegmentLabel, pData(datos)$trial_setting, pData(datos)$tumor_setting, sep = "_")
replicas_por_grupo <- pData(datos) %>%
as.data.frame() %>%
count(Grupo_gmx, name = "n_AOIs") %>%
arrange(n_AOIs)
replicas_por_grupo %>%
ggplot(aes(x = reorder(Grupo_gmx, n_AOIs), y = n_AOIs)) +
geom_col(fill = "steelblue") +
coord_flip() +
geom_hline(yintercept = 3, linetype = "dashed", color = "red") +
theme_bw() +
labs(
x = "Anotación biológica",
y = "Número de AOIs",
title = "Número de réplicas (AOIs) por anotación"
)
Finalmente, se obtienen 481 AOI y 81 variables para empezar el análisis.
En este punto del pipeline se van a utilizar mayoritariamente
funciones de StandR. Por este motivo es necesario cambiar la clase de
datos de un objeto NanostringGeoMxSet a un
objeto SpatialExperiment para poder aplicarle las funciones
pertinentes.
Para esta conversión, lo único que hay que tener en cuenta es que
SpatialExperiment se construye en “Target Level” y no en
“Probe level” que es como lo construye NanostringGeoMxSet.
Lo que significa que en la matriz de conteos aparecen genes únicos y no
sondas que pueden pertenecer al mismo gen. Se realiza la conversión de
los counts antes de pasarlo a un objeto
SpatialExperiment.
La clase SpatialExperiment tiene sus propias funciones
para acceder a cada parte del objeto, en la tabla siguiente se resume la
conversión de las funciones más destacadas.
- Conteos:
NanostringGeoMxSet->exprs(),SpatialExperiment->assays(). - Anotaciones AOI:
NanostringGeoMxSet->pData(),SpatialExperiment->colData(). - Anotaciones genes:
NanostringGeoMxSet->fData(),SpatialExperiment->rowData().
datos_probes <- datos
datos <- aggregateCounts(datos_probes)
datos_spe <- SpatialExperiment(
assays = list(counts = exprs(datos)),
rowData = fData(datos),
colData = pData(datos),
)Se añade logCPM(logCountsPerMilion) como
logcountsen datos, una métrica calculada a
partir de counts necesaria para pasos posteriores.
Es una métrica básica para normalizar por la profundidad de lectura y transformar con logaritmos para estabilizar la varianza.
logcpm <- cpm(assay(datos_spe, "counts"), log = TRUE, prior.count = 1)
assay(datos_spe, "logcounts") <- logcpm
assayNames(datos_spe)## [1] "counts" "logcounts"2.3.1 Diagrama de Sankey
Para la comprensión de la relación entre la matriz de conteos (counts por ROI) con cada muestra, se realiza un diagrama de Sankey. Este realiza la asociación entre la muestra, los pacientes, la clase (tratamiento o control), y el tipo celular (células tumorales (SOX2+) o TAM (IBA1+)).
Ejemplo 1 (GeomxTools):
library(dplyr)
library(ggforce)
library(networkD3)
library(htmlwidgets)
sankeyCols <- c("source", "target", "value")
link1 <- count(metadata, SlideName, trial_setting)
link2 <- count(metadata, trial_setting, tumor_setting)
link3 <- count(metadata, tumor_setting, SegmentLabel)
colnames(link1) <- sankeyCols
colnames(link2) <- sankeyCols
colnames(link3) <- sankeyCols
links <- rbind(link1,link2, link3)
nodes <- unique(data.frame(name=c(links$source, links$target)))
# sankeyNetwork is 0 based, not 1 based
links$source <- as.integer(match(links$source,nodes$name)-1)
links$target <- as.integer(match(links$target,nodes$name)-1)
sankey<- sankeyNetwork(Links = links, Nodes = nodes, Source = "source",
Target = "target", Value = "value", NodeID = "name",
units = "AOI", fontSize = 12, nodeWidth = 30)
sankeyEjemplo 2 (StandR):
plotSampleInfo(datos_spe, column2plot = c("SlideName", "trial_setting", "tumor_setting","SegmentLabel"))
3. Control de calidad (QC) y filtrado
El control de calidad en StandR se realiza a dos niveles: a nivel de gen y a nivel de ROI.
Antes de empezar con el QC se debe comprobar la columna
QCFlags que se encuentra en las anotaciones de los AOI
colData(datos_spe). Esta variable indica las AOI de mala
calidad según su paso preliminar del control de calidad y deben ser
eliminadas.
Si los datos superan correctamente su QC habrá valores NA o celdas vacías en esta columna.
## [1] TRUE## [1] 342## [1] 139qc_table <- as.data.frame(table(colData(datos_spe)$QCFlags))
names(qc_table) <- c("QC Flag", "Frecuencia")
kable(qc_table, caption = "Distribución de las banderas de control de calidad (QC)")| QC Flag | Frecuencia |
|---|---|
| Low Negative Probe Count for Probe Kit Human NGS Whole Transcriptome Atlas RNA_1.0 | 32 |
| Low Negative Probe Count for Probe Kit Human NGS Whole Transcriptome Atlas RNA_1.0,Low Nuclei Count | 2 |
| Low Negative Probe Count for Probe Kit Human NGS Whole Transcriptome Atlas RNA_1.0,Low Surface Area | 47 |
| Low Negative Probe Count for Probe Kit Human NGS Whole Transcriptome Atlas RNA_1.0,Low Surface Area,Low Nuclei Count | 18 |
| Low Nuclei Count | 3 |
| Low Raw Reads,Low Sequencing Saturation,Low Percent Aligned Reads,Low Percent Stitched Reads,Low Negative Probe Count for Probe Kit Human NGS Whole Transcriptome Atlas RNA_1.0,Low Percent Trimmed Reads | 0 |
| Low Surface Area | 33 |
| Low Surface Area,Low Nuclei Count | 4 |
Se obtienen 342 valores NA, que han superado el estándar de calidad.
Se obtienen 139 con avisos.
La tabla indica qué clases de QCFlags existen y cuantas
hay de cada clase.
Se filtran las AOIs con QCFlags que tengan valores
diferentes a NA.
## [1] 18677 342Se empieza el QC con 18677 genes y 342 AOI.
3.1 Gen QC
El control de calidad a nivel de gen tiene como objetivo eliminar aquellos genes que no se expresan en más del 90% de los AOI y detectar genes con expresión muy baja, menos de 5 counts. Es una forma de revisar la calidad de cada fila de la matriz de conteos.
Se calcula el umbral mínimo de expresión, expresado en logCPM, que
refleja el mínimo de counts esperado (5 en este caso) teniendo en cuenta
la profundidad de las librerías. Luego, se filtran los genes que no
alcanzan el min_count en más del 90% de las AOIs.
La función utilizada es addPerROI con los parámetros
rm_genes = TRUE (remove_genes), min_count = 5
y sample_fraction = 0.9.
datos_pre_QC <- datos_spe
datos_spe <- addPerROIQC(datos_pre_QC, rm_genes = TRUE)
names(metadata(datos_spe))## [1] "lcpm_threshold" "genes_rm_rawCount" "genes_rm_logCPM"Se generan tres archivos en el objeto datos_spe a los
cuales se accede mediante la función metadata():
lcpm_thresholdcontiene el umbral.genes_rm_rawCountcontiene los conteos crudos de los genes eliminados.genes_rm_logCPMcontiene los logcounts de los genes eliminados.
## [1] "NICN1" "DNAI3"Se eliminan 2 genes que son NICN1, DNAI3.
Con la función plotGeneQC() se evaluan las expresiones
en logCPM de los genes que fueron eliminados en las distintas AOIs.
La función también incluye un histograma que muestra la distribución del porcentaje de genes no expresados en todas las AOIs. Es decir, se calcula para cada AOI el porcentaje de genes que no se expresan (counts=0) y se hace un histograma de estos porcentajes.
Se obtienen dos gráficos, el primero muestra los genes eliminados, donde cada punto es la expresión de ese gen en un AOI (n=342). La mayoría de los puntos están por debajo del umbral de expresión (línea roja), lo que contribuye a que el gen sea eliminado. No se observa ninguna diferencia aparente entre los segmentos Iba1 y Sox2.
El segundo muestra para todos los AOI el porcentaje de genes que no tienen expresión. Se observa como la mayoría de AOI se acumulan mayoritariamente entre el 0-0,2% indicando que la mayoría de genes se expresan en todos los AOIs.
3.2 ROI QC
El control de calidad a nivel de ROI, busca identificar ROIs con bajo tamaño de biblioteca y/o bajo recuento de células, ya que se consideran muestras de baja calidad por insuficiente profundidad de secuenciación o falta de ARN en la región seleccionada. Es una forma de analizar la calidad de las columnas de la matriz de conteos. Como en este experimento se han recogido los conteos por AOI, se realiza el control de calidad de las AOI.
Se evalúa el gráfico de distribución del tamaño de la librería frente al número de núcleos. Al observar el diagrama de dispersión, se espera que los tamaños de librería estén correlacionados de forma positiva con el número de células (es decir, con el recuento de núcleos).
Para filtrar muestras de baja calidad, se define un umbral mínimo de recuento de células, para este experimento 150 células. Además, se puede investigar si los AOI de baja calidad provienen de alguna muestra específica, estratificando los puntos por nombre de lámina (slideName).
## `geom_smooth()` using formula = 'y ~ x'
No resulta inesperado que existan algunos AOIs con un tamaño de librería relativamente bajo pero con un número razonable de células (recuento de núcleos).
En este dataset, la relación entre tamaño de biblioteca y recuento de células es relativamente estable, sin anomalías evidentes, exceptuando una pequeña desviación al inicio de la recta. Esto indica que los valores observados son más altos de lo esperado en librerías con bajo número de recuento de células.
La línea roja marca el umbral de recuento mínimo de células = 150. No se observa ninguna asociación entre AOIs de baja calidad y una muestra en específico.
Eliminamos las AOI que no superan el número mínimo de células.
## qc
## FALSE TRUE
## 20 322## [1] 18675 322Se han eliminado 20 AOI.
Se pueden cambiar los ejes x o y de la función
plotROIQCpara cualquier otro estadístico que se desee
evaluar, usando los parámetros x_axis o y_axis.
Por ejemplo, se puede graficar el tamaño del área frente al tamaño de la biblioteca.
plotROIQC(datos_spe, x_axis = "AOISurfaceArea", x_lab = "AreaSize", y_axis = "lib_size", y_lab = "Library size", col = SlideName)## `geom_smooth()` using formula = 'y ~ x'
Se observa la presencia de un outlier, un punto que se encuentra en SaraWTATMA5 con un tamaño de librería entre 100.000-125.000 mm2. Se procede a eliminarlo.
outlier <- colData(datos_spe)$AOISurfaceArea >100000
datos_spe_outlier <- datos_spe
datos_spe <- datos_spe_outlier[, !outlier]
dim(datos_spe)## [1] 18675 321Se realiza el mismo gráfico sin el outlier.
plotROIQC(datos_spe, x_axis = "AOISurfaceArea", x_lab = "AreaSize", y_axis = "lib_size", y_lab = "Library size", col = SlideName)## `geom_smooth()` using formula = 'y ~ x'
La eliminación del outlier proporciona una relación entre el tamaño de librería y el tamaño del área de las AOI mucho más estable. Se observa claramente una relación positiva entre estos parámetros, indicando que a mayor tamaño de área de las AOI, mayor es la librería.
3.3 Expresión Logarítmica Relativa (RLE)
Tras el filtrado, se emplea la función plotRLExpr para
visualizar la Expresión Logarítmica Relativa (RLE) de los datos, con el
objetivo de identificar posibles variaciones técnicas presentes en el
conjunto de datos.
Para cada gen se calcula su mediana de expresión en todas las AOI. Luego se resta este valor a la expresión del gen de cada muestra.
Para cada gen:
RLE = log(expr_AOI) − mediana(log(expr_todas_las_AOI))
Es decir, el RLE indica cuánto se desvía la expresión de un gen en una AOI respecto al valor típico de ese gen. Para ello, se analiza la distancia relativa entre la mediana del RLE de cada AOI (representada por el punto en el boxplot) y el valor \(0\).
Por defecto, se representa el RLE de los datos crudos, donde se espera que la mayor parte de la variación observada esté causada por diferencias en el tamaño de la librería. Por lo tanto, el hecho de normalizar los datos solucionaría este problema.
Se puede observar que la mayoría de las AOI están centradas en la línea roja (\(y=0\)), aunque también se observan algunas desplazadas hacia arriba o con bigotes muy largos.
Si se realiza el RLE con los conteos logcounts, se encuentran en
assay=2, (recordatorio: logcounts es la transformación de
counts a logCPM), se observa como las variaciones técnicas debidas al
tamaño de librería son prácticamente eliminadas.
Existe la posibilidad de estratificar por anotaciones del metadata,
como por ejemplo ordannots = SlideName o
SegmentLabel, para identificar posibles factores que
contribuyen a las variaciones técnicas observadas.
No se encuentran grandes diferencias aparentemente entre los dos
grupos de SegmentLabel ni entre las muestras
SlideName.
3.4 Reducción de dimensionalidad
3.4.1 Análisis de Componentes Principales (PCA)
El PCA es un método de reducción de la dimensionalidad utilizado para transformar conjuntos de datos de alta dimensionalidad a dos dimensiones, sin perder la información más relevante. Ayuda a identificar agrupamientos (clústeres) de muestras con perfiles biológicos parecidos y a visualizar posibles variaciones técnicas.
Se realiza el PCA estratificando por SlideName para
evaluar si hay efecto de lote.
## Cargando paquete requerido: scuttle##
## Adjuntando el paquete: 'scater'## The following object is masked from 'package:limma':
##
## plotMDS## The following object is masked from 'package:standR':
##
## plotMDSset.seed(100)
datos_spe <- scater::runPCA(datos_spe)
pca_results <- reducedDim(datos_spe, "PCA")
drawPCA(datos_spe, precomputed = pca_results, col = SlideName)
El primer componente (horizontal) explica el 38.07% de la variación de los datos, mientras que el segundo (vertical) el 8.57%. Esta reducción a dos dimensiones explica el 46.64% de la variabilidad de los datos.
Por lo que se puede observar, no hay una agrupación clara de los
datos según el nombre de muestra, ya que no forman grupos aislados por
color. La mayor diferencia (PC1) parece ser de origen biológico o de
otro factor técnico que no es el SlideName, ya que se
intuyen dos grupos en el eje horizontal.
Se estratifica por el tipo de marcador (Iba1 o Sox2) y la clase (Nivolumab o control).
El análisis de componentes principales (PCA) revela que la principal fuente de variación en los datos (PC1: 38.07%) está impulsada por el tipo de segmento (Iba1 vs Sox2). No se observa una separación clara basada en el tratamiento (Nivolumab vs Control) a nivel global, ni un efecto de lote significativo por portaobjetos (SlideName), lo cual valida la consistencia biológica de las muestras.
Otras funciones del paquete standR para la visualización de los PCAs.
## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
3.4.2 Multidimensional Scaling (MDS)
Otra forma de reducir dimensionalidad es mediante el MDS.
standR::plotMDS(datos_spe,
dims = c(1, 2),
precomputed = NULL,
textScale = 1,
assay = 1,
color = SegmentLabel)
Se explica una mayor variabilidad de los datos del 60.56% en estas dos dimensiones, aunque la separación entre los tipos de segmento Iba1 y Sox2 es menos clara que cuando se usa PCA.
3.4.3 Uniform Manifold Approximation and Projection (UMAP)
Otro método de reducción de la dimensionalidad es mediante UMAP.
set.seed(100)
datos_spe <- scater::runUMAP(datos_spe, dimred = "PCA")
plotDR(datos_spe, dimred = "UMAP", col = SegmentLabel)
Se corrobora que la presencia de variabilidad aportada por el efecto de lote es mínima y que la separación entre grupos es debida a las dos poblaciones Iba1 y Sox2.
4. Normalización
Si se identifican variaciones técnicas en los pasos de QC, como las observadas en el RLE de los datos crudos, antes de realizar cualquier análisis posterior es necesario normalizar los datos para corregir o minimizar dichas variaciones.
El paquete standR ofrece varias opciones de normalización: TMM, RPKM, TPM, CPM, upperquartile y sizefactor.
Se realiza la normalización por el método TMM mediante la función
geomxNorm(). Esta función normaliza los counts
de datos_spe, y almacena la matriz normalizada en un nuevo
assay bajo el nombre de logcounts.
De esta forma, se obtiene en datos_spe y
datos_tmm, la misma matriz counts y matrices
distintas en logcounts, ya que en datos_spe,
logcounts = logCPM, en cambio, en datos_tmm,
logcounts= datos normalizados por TMM.
## [1] "counts" "logcounts"Para evaluar si la normalización eliminó las variaciones no deseadas, se utilizan los gráficos RLE y PCA.
library(patchwork)
p1 <- plotRLExpr(datos_tmm, ordannots = "SlideName", assay = 1, color = SlideName) + ggtitle ("RAW")
p2 <- plotRLExpr(datos_tmm, ordannots = "SlideName", assay =2, color = SlideName) + ggtitle("TMM")
p1 
Se observa como la normalización aporta una gran mejora acercando las medianas de las AOI a 0, sugiriendo que la eliminación de las variaciones técnicas.
set.seed(100)
datos_tmm <- scater::runPCA(datos_tmm, exprs_values = "logcounts")
pca_results_tmm <- reducedDim(datos_tmm, "PCA")
plotPairPCA(datos_tmm, precomputed = pca_results_tmm, color = SegmentLabel)
El PC1 (37.09%) sigue siendo el eje principal de variación y separa
casi perfectamente los dos grupos biológicos (Iba1 y Sox2). Aunque el
porcentaje de varianza del PC1 baja ligeramente (de 38.07% a 37.09%), se
ha eliminado parte del “ruido técnico” que previamente inflaba
artificialmente la varianza. Si se compara la distribución por
SlideName, los colores de los diferentes portaobjetos están
completamente mezclados entre sí. Esto confirma que la normalización TMM
ha eliminado con éxito el efecto de lote.
5. Corrección del efecto lote
Cada lámina generalmente solo puede contener unos pocos segmentos de tejido, por lo que es común que los datos de GeoMx DSP estén afectados por el efecto de lote introducido por las diferentes láminas.
En los pasos previos a la corrección del efecto lote, se ha encontrado una presencia de efecto de lote mínima en este dataset. Igualmente, es recomendable realizar la corrección del efecto lote, para augmentar el poder estadístico, para minimizar los falsos positivos (que un gen aparezca como significativo cuando en verdad no lo es ) y porque se considera una buena praxis demostrando robusted en el análisis.
El paquete standR ofrece dos enfoques para corregir efectos de lote: RUV4 y Limma.
5.1 RUV4
La función findNCGs permite identificar los NCGs (Negative Control Genes) de los datos. En este caso, el efecto de lote se debe principalmente a las láminas, así que se busca identificar NCGs a través de todas las láminas.
Para ello, el parámetro batch_name se establece en “SlideName” y se seleccionan los 300 genes menos variables de las distintas láminas como NCGs (ordenados por coeficiente de variación). Se hace esta selección de genes estables para tener una referencia neutra, así cualquier variación observada en los demás genes reflejará principalmenteuna duda diferencias biológicas y no efectos de lote.
Estos genes se almacenan en el metadata(datos_spe) como
una columna nueva bajo el nombre “NCGs”.
La función findNCG() utiliza por defecto el
assay = 2 = logcounts de los datos introducidos.
## New names:
## • `cv` -> `cv...1`
## • `cv` -> `cv...2`
## • `cv` -> `cv...3`
## • `cv` -> `cv...4`
## • `cv` -> `cv...5`
## • `cv` -> `cv...6`## [1] "lcpm_threshold" "genes_rm_rawCount" "genes_rm_logCPM"
## [4] "NCGs"Una de las principales dudas a la hora de aplicar la corrección del
efecto lote del pipeline StandR, es el por qué se aplica a los
logcounts de datos_spe y no a los
logcounts de datos_tmm.
Por este motivo se realiza también la corrección del efecto lote en los datos normalizados.
## New names:
## • `cv` -> `cv...1`
## • `cv` -> `cv...2`
## • `cv` -> `cv...3`
## • `cv` -> `cv...4`
## • `cv` -> `cv...5`
## • `cv` -> `cv...6`## [1] "lcpm_threshold" "genes_rm_rawCount" "genes_rm_logCPM"
## [4] "norm.factor" "norm.method" "NCGs"Para la corrección con RUV4, la función requiere 3 parámetros adicionales además del objeto de entrada:
factors: el factor de interés, es decir, la variación biológica que se desea conservar.NCGs: la lista de genes de control negativo identificados con la función findNCGs.k: el número de factores no deseados a utilizar. Según la documentación de RUV, se recomienda usar el \(k\) más pequeño posible que elimine la variación técnica observada.
La función geomBatchCorrection() también utiliza por
defecto el assay = 2 = logcounts
Se inspecciona sobre qué \(k\) es la mejor para este análisis. El valor óptimo de \(k\) será el valor más pequeño que logre producir una separación de la biología principal de interés del experimento en un gráfico de PCA.
for(i in seq(5)){
datos_ruv <- geomxBatchCorrection(datos_spe, factors = "SegmentLabel",
NCGs = metadata(datos_spe)$NCGs, k = i)
print(plotPairPCA(datos_ruv, assay = 2, n_dimension = 4, color = SegmentLabel, title = paste0("k = ", i)))
}
Se elije \(k=2\) como mejor resultado.
datos_ruv <- geomxBatchCorrection(datos_spe, factors = "SegmentLabel",
NCGs = metadata(datos_spe)$NCGs, k = 2)
set.seed(100)
datos_ruv <- scater::runPCA(datos_ruv)
pca_results_ruv <- reducedDim(datos_ruv, "PCA")
plotPairPCA(datos_ruv, precomputed = pca_results_ruv, color = SegmentLabel, title = "RUV4, k = 2", n_dimension = 4)
plotPairPCA(datos_ruv, precomputed = pca_results_ruv, color = SlideName, title = "RUV4, k = 2", n_dimension = 4)
Se observa que se ha producido la eliminación del efecto de lote,
porque no se puede intuir ninguna agrupación en el PCA estratificado por
SlideName.
Se repite la misma labor para los datos normalizados.
for(i in seq(5)){
datos_ruv_norm <- geomxBatchCorrection(datos_tmm, factors = "SegmentLabel",
NCGs = metadata(datos_tmm)$NCGs, k = i)
print(plotPairPCA(datos_ruv_norm, assay = 2, n_dimension = 4, color = SegmentLabel, title = paste0("k = ", i)))
}
Se elije \(k=2\) como valor óptimo.
datos_ruv_norm <- geomxBatchCorrection(datos_tmm, factors = "SegmentLabel",
NCGs = metadata(datos_tmm)$NCGs, k = 2)
set.seed(100)
datos_ruv_norm <- scater::runPCA(datos_ruv_norm)
pca_results_ruv_norm <- reducedDim(datos_ruv_norm, "PCA")
plotPairPCA(datos_ruv_norm, precomputed = pca_results_ruv_norm, color = SegmentLabel, title = "RUV4, k = 2", n_dimension = 4)
plotPairPCA(datos_ruv_norm, precomputed = pca_results_ruv_norm, color = SlideName, title = "RUV4, k = 2", n_dimension = 4)
Los PCAs de los datos normalizados muestran una reducción del ruido
técnico, logrando una integración de los diferentes portaobjetos
SlideName. Se obtiene una subida de la explicación del PC2
a 11.19%, sugiriendo que existen diferencias secundarias dentro de cada
tipo celular.
Comparando los PCAs de las primeras dos dimensiones de los
logcounts de los datos crudos y normalizados, se percibe
una corrección más limpia del efecto de lote en los datos crudos
datos_spe, ya que se observa una mejor dispersión de los
datos según su slideName.
5.2 Limma
La segunda opción en el pipeline StandR para la eliminación del ruido de fondo es el paquete Limma. Esta función requiere 2 parámetros adicionales además del objeto de entrada.
batch: un vector que indica la información de lote de todas las muestras.design: una matriz de diseño generada a partir demodel.matrix.
datos_lrb <- geomxBatchCorrection(datos_spe,
batch = colData(datos_spe)$SlideName, method = "Limma",
design = model.matrix(~SegmentLabel, data = colData(datos_spe)))Nuevamente se utilizan gráficos de QC, como el PCA, para inspeccionar y evaluar la efectividad del proceso de corrección de efectos de lote aplicado.
## Warning in (function (mapping = NULL, data = NULL, stat = "identity", position
## = "identity", : Ignoring unknown parameters: `title`
Se visualiza la corrección del efecto lote, ya que se observan los puntos de las distintas muestras aleatoriamente repartidos y sin ningún patrón de agrupación.
Se realiza, al igual que con RUV4, el análisis con los datos normalizados.
datos_lrb_norm <- geomxBatchCorrection(datos_tmm,
batch = colData(datos_tmm)$SlideName, method = "Limma", n_assay = 2,
design = model.matrix(~SegmentLabel, data = colData(datos_tmm)))
plotPairPCA(datos_lrb_norm, assay = 2, color = SegmentLabel, title = "Limma remove Batch Norm TMM Data")
plotPairPCA(datos_lrb_norm, assay = 2, color = SlideName, title = "Limma remove Batch Norm TMM Data")
## Warning in (function (mapping = NULL, data = NULL, stat = "identity", position
## = "identity", : Ignoring unknown parameters: `title`
Se observa en estos últimos gráficos, una separación menos nítida
estratificando por SegmentLabelen comparación con los datos
no normalizados, cuando se aplica la corrección del efecto lote con
Limma. Esto no significa la obtención de peores resultados, sino que
indica que la separación observada en los logcountsde
datos_spe estaba parcialmente inflada. Tras la corrección
se obtiene una visión más fiel y conservadora de la verdadera
variabilidad biológica entre los segmentos.
Se recomienda utilizar estadísticas resumidas para evaluar la efectividad de la corrección por lote y elegir qué método (RUV4 o Limma) ha sido más efectivo para este dataset.
Las 6 estadísticas resumidas incluidas en el paquete son conocidas como Métricas de Evaluación de Clústeres:
Adjusted Rand Index (ARI): Mide la concordancia entre los clústeres del PCA y las etiquetas reales (SegmentLabel), ajustando el valor por el azar (donde 1 es una coincidencia perfecta).
Coeficiente de similitud de Jaccard: Evalúa qué tan similares son los conjuntos de muestras en cada clúster, calculando la proporción de coincidencias sobre el total de elementos comparados.
Coeficiente de Silhouette: Mide qué tan cerca está cada muestra de su propio grupo en comparación con otros grupos; valores altos indican una asignación de clúster clara y robusta.
Coeficiente Chi-cuadrado: Evalúa la independencia estadística entre los clústeres obtenidos y las variables de lote o segmento, detectando si existe un sesgo sistemático.
Distancia de Mirkin: Cuantifica la diferencia o “distancia” entre dos clasificaciones; a menor valor, mayor es la similitud entre el agrupamiento del PCA y la biología esperada.
Coeficiente de solapamiento (Overlap Coefficient): Mide el grado de “intersección” entre los grupos, siendo muy útil para comprobar si las muestras de diferentes lotes se han integrado correctamente en un mismo espacio.
Esta evaluación se realiza con la función
plotClusterEvalStats().
Se presentan los resultados de las estadísticas resumidas para los dos métodos de corrección del efecto de lote utilizados en este pipeline (RUV4 y Limma), así como para los datos crudos y normalizados con TMM.
En orden siguiente: datos_spe, datos_ruv, datos_lrb, datos_tmm, datos_ruv_norm, datos_lrb_norm.
La puntuación de cada método se muestra en un gráfico de barras con las seis estadísticas resumidas, organizadas en dos secciones: biología y lote.
Para los factores biológicos de interés definidos en la corrección por lote, es mejor un puntaje más alto en la corrección.
Para los efectos de lote, un puntaje más bajo es preferible en la corrección.
datos_list <- list(datos_spe, datos_ruv, datos_lrb, datos_tmm, datos_ruv_norm, datos_lrb_norm)
plotClusterEvalStats(spe_list = datos_list,
bio_feature_name = "SegmentLabel",
batch_feature_name = "SlideName",
data_names = c("Raw","RUV4","Limma", "Raw TMM", "RUV4 TMM", "Limma TMM"))
En este conjunto de datos, se identifica RUV4 en los datos sin normalizar, como la herramienta más potente para preservar el carácter biológico de las muestras y Limma en los datos normalizados, como la mejor herramienta para minimizar el ruido técnico.
Para finalizar, se pueden visualizar mediante gráficos RLE, los
logcounts de datos_ruv/lrb y de
datos_ruv/lrb_norm para determinar qué corrección por lote
funciona mejor en este conjunto de datos.
Se observa un efecto similar ente ambos métodos, así que se puede usar indistintamente cualquiera de los dos para este conjunto de datos.
Al comparar los logcounts con los datos normalizados por TMM, se intuye una mayor centralidad de las cajas en los resultados normalizados. Esto nos indica que el hecho de normalizar comporta una mejoría a la hora de poder comparar muestras que se encuentran a diferente escala.
En este punto, elegir la metodología adecuada para la corrección, va
a ser influido por el objetivo del experimento, diferenciando en si la
prioridad es la precisión biológica (datos_ruv), o la
integración perfecta de las muestras (datos_lrb_norm).
Se elije el método RUV4 (datos_ruv) para continuar con
los análisis, priorizando la precisión biológica.
6. Análisis posteriores (downstream)
Para los análisis posteriores, como los de expresión diferencial (DE), standR no proporciona funciones específicas.
Se recomienda integrar el flujo de trabajo con pipelines establecidos, como edgeR, limma-voom o DESeq2, que utilizan modelos lineales y aprovechan la información de todos los genes, siendo más adecuados para datasets complejos con múltiples factores experimentales. No se recomienda usar análisis simples como t-test para analizar datos de GeoMx DSP.
En este pipeline, se desarrollará el análisis DE usando el pipeline limma-voom.
Se ha demostrado previamente que, para este dataset, el uso de RUV
con \(k=2\) (datos_ruv)
constituye la estrategia más equilibrada; corrigeeficazmente el efecto
de lote y otras variaciones técnicas no deseadas sin comprometer la
resolución de la estructura biológica de los datos.
Sin embargo, los conteos normalizados no se deben usar directamente
en modelos lineales. Para el modelado lineal es mejor incluir la matriz
de pesos generada por la función geomxBatchCorrection()
como covariable. Esta matriz de pesos se encuentra en el
colData() del objeto.
pesos_ruv <- colData(datos_ruv)[,seq(ncol(colData(datos_ruv))-1, ncol(colData(datos_ruv)))]
head(pesos_ruv)## DataFrame with 6 rows and 2 columns
## ruv_W1 ruv_W2
## <numeric> <numeric>
## DSP-1001660012264-D-A04.dcc -0.0183241 0.13944120
## DSP-1001660012264-D-A06.dcc -0.0147220 -0.06229172
## DSP-1001660012264-D-A07.dcc -0.0428157 -0.00670244
## DSP-1001660012264-D-A09.dcc -0.0129535 -0.06966592
## DSP-1001660012264-D-A10.dcc -0.0204593 -0.06208149
## DSP-1001660012264-D-B06.dcc 0.0421647 -0.054052886.1 Establecimiento de una matriz de diseño y contrastes
El objetivo del análisis posterior es comprender el carácter biológico de los datos que se están analizando. Para ello se define como objetivo evaluar si existe una diferencia significativa entre cómo se expresan las células cancerígenas (Sox2) y los macrófagos asociados a tumores (Iba1) en los pacientes control y tratados con neoadyuvante.
En teoría, si el tratamiento funciona y hay una activación de la resupuesta inmunitária, esto debería reflejarse en cambios de expresión génica en las distintas zonas celulares analizadas. Un ejemplo de este cambio de perfil de expresión podría ser:
En células tumorales: aumento de genes de estrés, apoptosis, interferones.
En TAMs: activación de pro-inflamatoria (M1), desactivación de inmunosupresores, secreción de citoquinas y quimiocinas.
También se desea analizar si existen diferencias significativas entre la eficacia del tratamiento en pacientes con tumores recurrentes o primarios. Los tumores primarios son potencialmente más susceptibles al tratamiento, mientras que los tumores recurrentes tienen un microambiente preparado para resistirlo. Por lo tanto, pueden existir diferencias en la activación de genes en los pacientes tratados con Nivolumab según si el tumor es primario o recurrente.
6.2.1 Análisis de Expresión Diferencial (DE) Iba1 vs Sox2
Para iniciar el análisis se realiza un enfoque en la expresión diferencial inicial entre los segmentos Iba1 y Sox2, para observar que los perfiles transcriptómicos sean distintos, lo cual es esperado porque son tipos celulares diferentes.
Se construye una matriz de diseño con la información de los subtipos de tejido.
Se añaden las matrices W resultantes de RUV4 a la matriz de diseño como covariables para usar los datos corregidos por efecto de lote.
dge <- SE2DGEList(datos_ruv)
design1 <- model.matrix(~0 + SegmentLabel + ruv_W1 + ruv_W2, data = colData(datos_ruv))
colnames(design1)## [1] "SegmentLabelIba1" "SegmentLabelSox2" "ruv_W1" "ruv_W2"Se redefinen los nombres.
## [1] "Iba1" "Sox2" "ruv_W1" "ruv_W2"Se construye la matriz de contraste para el análisis de Expresion Diferencial (DE) entre los segmentos.
Se recomienda filtrar los genes con baja cobertura en el dataset para obtener una relación media-varianza más precisa y reducir el número de pruebas estadísticas.
En este caso, se utiliza la función filterByExpr del
paquete edgeR para filtrar los genes según la matriz de
diseño, conservando la mayor cantidad posible de genes con conteos
razonables.
## keep
## FALSE TRUE
## 2 18673## [1] "NPBWR2" "SPATC1"Aparecen solo 2 genes de baja cobertura, que son NPBWR2, VWC2. Se eliminan.
6.2.2 Comprovación del Coeficiente de Variación Biológico (BCV)
El Coeficiente de Variación Biológico (BCV) es un indicador de la variabilidad biológica entre replicados o muestras biológicas. Refleja cómo cambia la abundancia real (desconocida) de un gen entre muestras replicadas de ARN.
En el gráfico de BCV de un dataset de GeoMx hay tres aspectos principales a observar:
Tendencia de dispersión: Se espera que la dispersión se aplane en los genes con mayor conteo.
Tendencia común: Debe ser relativamente pequeña. En RNA-seq, se espera entre 0.2 y 0.4 en muestras humanas y entre 0.05 y 0.2 en ratón y líneas celulares. En GeoMx, al muestrear segmentos de tejido humano, se espera que sea menor que en RNA-seq humano.
Genes con bajo conteo y alto BCV: Hay que tener cuidado si aparecen genes con alto BCV y bajo conteo. A la vez es una buena praxis revisar los genes con alto BCV con especialistas, ya que probablemente se identifiquen como genes diferencialmente expresados (DE) y puedan estar influyendo en la variación observada en los gráficos PCA.
Se genera el gráfico de BCV para la comparación entre los segmentos Iba1 y Sox2.
dge_segment <- estimateDisp(dge_filtrado, design = design1, robust = TRUE)
plotBCV(dge_segment, legend.position = "topleft", ylim = c(0, 1.3))## Warning in plot.window(...): "legend.position" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "legend.position" es un parámetro gráfico
## inválido## Warning in axis(side = side, at = at, labels = labels, ...): "legend.position"
## es un parámetro gráfico inválido
## Warning in axis(side = side, at = at, labels = labels, ...): "legend.position"
## es un parámetro gráfico inválido## Warning in box(...): "legend.position" es un parámetro gráfico inválido## Warning in title(...): "legend.position" es un parámetro gráfico inválidobcv_df <- data.frame(
'BCV' = sqrt(dge_segment$tagwise.dispersion),
'AveLogCPM' = dge_segment$AveLogCPM,
'gene_id' = rownames(dge_segment)
)
highbcv <- bcv_df$BCV > 0.8
highbcv_df <- bcv_df[highbcv, ]
points(highbcv_df$AveLogCPM, highbcv_df$BCV, col = "red")
text(highbcv_df$AveLogCPM, highbcv_df$BCV, labels = highbcv_df$gene_id, pos = 4)
La línea de tendencia de dispersión (azul) se acaba aplanando en los genes con mayores conteos. La línea de tendencia común (roja) está por debajo del rango humano de RNA-seq (0,2-0,4). No aparecen genes con alto BCV y bajo número de counts.
Los genes como GFAP, SPP1, CHI3L1 o PDPN tienen un coeficiente de variación biológico alto, lo que implica que su expresión varía mucho entre los segmentos celulares.
6.2.3 Análisis de Expresión Diferencial (DE) con el pipeline limma-voom
En el pipeline de análisis limma-voom, el modelado lineal se realiza sobre los valores log-CPM utilizando las funciones voom, lmFit, contrasts.fit y eBayes.
El método voom (de limma) estima cómo cambia la variabilidad de los genes dependiendo de su expresión. La línea roja debería suavizar la nube de puntos. Eso significa que voom estimó bien la relación media–varianza.
Cada punto negro es un gen.
En el eje X está el log2 del tamaño de los conteos+0,5.
En el eje Y está la desviación estándar estimada para cada gen.
La línea roja es la tendencia media-variancia ajustada.
Se obtiene una especie de U, donde a la izquierda se acumulan los genes poco expresados, en el centro hay una nube de puntos donde se encuentran los genes con expresión media (log2(count+0,5)=5/6) y a la derecha parecen genes muy expresados. Este gráfico sirve para poder evaluar si voom detecta una relación media-varianza clara. Si la línea roja sigue claramente la nube de puntos, voom está capturando bien la relación media–varianza y el modelo está bien ajustado. Con el resultado obtenido, se puede afirmar que el modelo estimó bien a relación media-varianza.
Se prosigue al ajuste del modelo con limma, donde
lmFit() ajusta el modelo, contraste.fit()
define comparaciones y eBayes aporta robustez en el cálculo
estadístico.
fit1 <- lmFit(v1)
fit_contrast1 <- contrasts.fit(fit1, contrasts = contr_matrix_segment)
efit1 <- eBayes(fit_contrast1, robust = TRUE)Se mide qué genes son estadísticamente significativos.
results_efit1<- decideTests(efit1, p.value = 0.05)
summary_efit1 <- summary(results_efit1)
summary_efit1## IvsS
## Down 5282
## NotSig 4976
## Up 8415Iba1 vs Sox2:
- Down: Genes con mayor expresión en Sox2 5282.
- NotSig: Genes con expresión similar en ambos segmentos 4976.
- Up: Genes con mayor expresión en Iba1 8415.
6.2.4 Visualización
Mediante TopTable() se obtienen los genes significativos
ordenados por pValue del análisis de Expresión
Diferencial.
DE_results_IvS <- topTable(efit1, coef = 1, sort.by = "P", n = Inf)
DE_genes_toptable_IvS <- topTable(efit1, coef = 1, sort.by = "P", n = Inf, p.value = 0.05)Se grafica con con un MA.
DE_results_IvS %>%
mutate(DE = ifelse(logFC > 0 & adj.P.Val <0.05, "UP",
ifelse(logFC <0 & adj.P.Val<0.05, "DOWN", "NOT DE"))) %>%
ggplot(aes(AveExpr, logFC, col = DE)) +
geom_point(shape = 1, size = 1) +
geom_text_repel(data = DE_genes_toptable_IvS %>%
mutate(DE = ifelse(logFC > 0 & adj.P.Val <0.05, "UP",
ifelse(logFC <0 & adj.P.Val<0.05, "DOWN", "NOT DE"))) %>%
rownames_to_column(), aes(label = rowname), max.overlaps = 20) +
theme_bw() +
xlab("Average log-expression") +
ylab("Log-fold-change") +
ggtitle("Iba1 vs. Sox2 (limma-voom)") +
scale_color_manual(values = c("blue","gray","red")) +
theme(text = element_text(size=15))## Warning: ggrepel: 13669 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
La matriz de contraste se calculó como
Contr.matrix= Iba1 - Sox2, así que se obtiene en rojo los
genes sobreexpresados en áreas Iba1, que corresponden a células
cancerosas y en azul los genes sobreexpresados en áreas Sox2, que
corresponden a macrófagos asociados a tumores. La zona gris son genes no
diferencialmente expresados entre ambas áreas.
Otra forma de acceder a la información es a través de una tabla interactiva usando el paquete DT.
updn_cols <- c(RColorBrewer::brewer.pal(6, 'Greens')[2], RColorBrewer::brewer.pal(6, 'Purples')[2])
DE_genes_toptable_IvS %>%
dplyr::select(c("logFC", "AveExpr", "P.Value", "adj.P.Val")) %>%
DT::datatable(caption = 'Iba1 zone vs. Sox2 zone in GlioBastomas (limma-voom)') %>%
DT::formatStyle('logFC',
valueColumns = 'logFC',
backgroundColor = DT::styleInterval(0, rev(updn_cols))) %>%
DT::formatSignif(1:4, digits = 4)## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.htmlSe deduce, tanto del MA plot como de la tabla, que genes como LAPTM5, ARGHGDIB, TYROBP, C1QC y SLCO2B1, son genes sobreexpresados en Iba1 en comparación con Sox2 y, si se evalúan sus funciones, son característicos de células mieloides, especialmente de macrófagos y microglía, lo cual encaja perfectamente con las células TAM.
6.3.1 Análisis de Expresión Diferencial (DE) para la variable Grupo: Iba1/Sox2 + tratamiento/control + Primary/Recurrence
Para poder hacer el análisis de los objetivos descritos, se genera una nueva matriz de diseño que contenga una variable compuesta por la información necesaria para luego poder realizar los análisis. Así que se construye una nueva variable compuesta por el SegmentLabel+ trial_setting + tumor_setting.
colData(datos_ruv)$grupo <- paste(colData(datos_ruv)$SegmentLabel, colData(datos_ruv)$trial_setting, colData(datos_ruv)$tumor_setting, sep = "_")Se construye una matriz de diseño con la variable grupo
que se compone de las 3 variables agrupadas, usando los datos donde se
ha eliminado el efecto de lote con el método Ruv4.
Se añaden las matrices W resultantes de RUV4 a la matriz de diseño como covariables.
## [1] "grupoIba1_Control_Primary" "grupoIba1_Control_Recurrence"
## [3] "grupoIba1_Nivolumab_Primary" "grupoIba1_Nivolumab_Recurrence"
## [5] "grupoSox2_Control_Primary" "grupoSox2_Control_Recurrence"
## [7] "grupoSox2_Nivolumab_Primary" "grupoSox2_Nivolumab_Recurrence"
## [9] "ruv_W1" "ruv_W2"Se redefinen los nombres.
## [1] "Iba1_Control_Primary" "Iba1_Control_Recurrence"
## [3] "Iba1_Nivolumab_Primary" "Iba1_Nivolumab_Recurrence"
## [5] "Sox2_Control_Primary" "Sox2_Control_Recurrence"
## [7] "Sox2_Nivolumab_Primary" "Sox2_Nivolumab_Recurrence"
## [9] "ruv_W1" "ruv_W2"Se obtienen 8 niveles, (2^3).
Se construye una matriz de contraste para el análisis de expresion diferencial (DE) entre tratamiento y control.
contr_matrix_tratamiento <- makeContrasts(
Iba1_Nivo_vs_Control_Rec = Iba1_Nivolumab_Recurrence - Iba1_Control_Recurrence,
SOX2_Nivo_vs_Control_Rec = Sox2_Nivolumab_Recurrence - Sox2_Control_Recurrence,
levels = colnames(design2)
)La expresión diferencial se evalúa utilizando estadísticos moderados empírico-bayesianos implementados en limma.
keep <- filterByExpr(dge, design2)
dge_filtrado2 <- dge[keep, ]
v2 <- voom(dge_filtrado2, design2, plot = TRUE)
fit2 <- lmFit(v2)
fit_contrast2 <- contrasts.fit(fit2, contrasts = contr_matrix_tratamiento)
efit2 <- eBayes(fit_contrast2, robust = TRUE)Se mide qué genes son estadísticamente significativos.
results_efit2<- decideTests(efit2, p.value = 0.05)
summary_efit2 <- summary(results_efit2)
summary_efit2## Iba1_Nivo_vs_Control_Rec SOX2_Nivo_vs_Control_Rec
## Down 266 355
## NotSig 17500 17321
## Up 907 997Se observa que la mayoría de genes forman parte de la clase No Significativos.
Se extraen los resultados de ambos contrastes en froma de tabla. Se genera un volcano plot, para una mejor observación de estos resultados.
resultado_Iba1 <- topTable(efit2, coef = "Iba1_Nivo_vs_Control_Rec", number = Inf, sort.by = "P")
resultado_Sox2 <- topTable(efit2, coef = "SOX2_Nivo_vs_Control_Rec", number = Inf, sort.by = "P")Se define el umbral de significación.
resultado_Iba1$significant <- with(resultado_Iba1, adj.P.Val < 0.05 & abs(logFC) > 1)
resultado_Sox2$significant <- with(resultado_Sox2, adj.P.Val < 0.05 & abs(logFC) > 1)Se crean los Volcano Plot.
library(ggplot2)
Volcano_Iba1 <- ggplot(resultado_Iba1, aes(x = logFC, y = -log10(adj.P.Val))) + geom_point(aes(color = significant), alpha = 0.6) +
scale_color_manual(values = c("grey", "red")) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
theme_minimal() +
labs(
title = "Iba1+ AOIs",
x = "log2 Fold Change",
y = "-log10(FDR)"
)
Volcano_Sox2 <- ggplot(resultado_Sox2, aes(x = logFC, y = -log10(adj.P.Val))) + geom_point(aes(color = significant), alpha = 0.6) +
scale_color_manual(values = c("grey", "blue")) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
theme_minimal() +
labs(title = "Sox2+ AOIs", x = "log2 Fold Change", y = "-log10(FDR)")
library(patchwork)
combined_volcano <- Volcano_Iba1 + Volcano_Sox2 + plot_layout(ncol = 2) + plot_annotation(title = "Volcano Plot de Expresión Diferencial", subtitle = "Nivolumab vs Control en pacietes Recurrentes")
combined_volcano
No se observan diferencias significativas entre la expresión génica de los pacientes recurrentes tratados con Nivolumab o pacientes recurrentes control, en ninguna de las dos zonas celulares analizadas. Esto nos indica que el tratamiento no está siendo efectivo, ya que no se está desarrollando una activación de la respuesta inmunitaria, que cambiaría el perfil de expresión genética de ambas zonas celulares.
A continuación se analizará si existe diferencia en la efectividad del tratamiento en personas con tumores primarios y recurrentes.
contr_matrix_tipotumor <- makeContrasts(
Iba1_Nivo_Pri_vs_Rec = Iba1_Nivolumab_Primary - Iba1_Nivolumab_Recurrence,
SOX2_Nivo_Pri_vs_Rec = Sox2_Nivolumab_Primary - Sox2_Nivolumab_Recurrence,
levels = colnames(design2)
)La expresión diferencial se evalua utilizando estadísticos moderados empírico-bayesianos implementados en limma.
Se mide qué genes son estadísticamente significativos.
fit_contrast3 <- contrasts.fit(fit2, contrasts = contr_matrix_tipotumor)
efit3 <- eBayes(fit_contrast3, robust = TRUE)
results_efit3<- decideTests(efit3, p.value = 0.05)
summary_efit3 <- summary(results_efit3)
summary_efit3## Iba1_Nivo_Pri_vs_Rec SOX2_Nivo_Pri_vs_Rec
## Down 197 52
## NotSig 18217 18563
## Up 259 58Se observa de nuevo como la mayoría de genes forman parte de la clase No Significativos.
Se extraen los resultados de ambos contrastes en froma de tabla.
Se define el umbral de significación.
Se genera un volcano plot, para una mejor observación de estos resultados.
resultado_Iba1T <- topTable(efit3, coef = "Iba1_Nivo_Pri_vs_Rec", number = Inf, sort.by = "P")
resultado_Sox2T <- topTable(efit3, coef = "SOX2_Nivo_Pri_vs_Rec", number = Inf, sort.by = "P")
resultado_Iba1T$significant <- with(resultado_Iba1T, adj.P.Val < 0.05 & abs(logFC) > 1)
resultado_Sox2T$significant <- with(resultado_Sox2T, adj.P.Val < 0.05 & abs(logFC) > 1)
Volcano_Iba1T <- ggplot(resultado_Iba1T, aes(x = logFC, y = -log10(adj.P.Val))) + geom_point(aes(color = significant), alpha = 0.6) +
scale_color_manual(values = c("grey", "red")) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
theme_minimal() +
labs(
title = "Iba1+ AOIs",
x = "log2 Fold Change",
y = "-log10(FDR)"
)
Volcano_Sox2T <- ggplot(resultado_Sox2T, aes(x = logFC, y = -log10(adj.P.Val))) + geom_point(aes(color = significant), alpha = 0.6) +
scale_color_manual(values = c("grey", "blue")) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
theme_minimal() +
labs(title = "Sox2+ AOIs", x = "log2 Fold Change", y = "-log10(FDR)")
library(patchwork)
combined_volcano <- Volcano_Iba1T + Volcano_Sox2T + plot_layout(ncol = 2) + plot_annotation(title = "Volcano Plot de Expresión Diferencial", subtitle = "Primary vs Recurrence en pacietes tratados con Nivolumab")
combined_volcano
Se observa que no hay diferencias significativas entre los pacientes con tumores primarios y recurrentes tratados con Nivolumab en ninguna de las dos zonas celulares. Esto demuestra que el tratamiento no es más efectivo en pacientes con tumores primarios, hecho que apoya de nuevo la idea de que el tratamiento no es efectivo.
Debido a los resultados encontrados hasta ahora, se decide realizar la comparación entre las zonas celulares Iba1 y Sox2. Para ello, se construye una nueva matriz que contrasta los grupos Iba1_Control_Primary vs Sox2_Control_Primary. Intentando representar condiciones “vírgenes” para analizar las diferencias entre las zonas celulares Iba1 y Sox2.
Se elige el método de corrección de efecto de lote Limma, a diferencia de Ruv4 implementado anteriormente, para familiarizarse también con este método en el pipeline.
colData(datos_lrb)$grupo <- paste(colData(datos_lrb)$SegmentLabel, colData(datos_lrb)$trial_setting, colData(datos_lrb)$tumor_setting, sep = "_")Se construye una matriz de diseño con la variable grupo
que se compone de las 3 variables agrupadas, usando los datos donde se
ha eliminado el efecto de lote con el método Limma.
En este caso no se añaden los pesos como covariables, sino que se debería añadir el factor que aporta variabilidad técnica como factor fijo en la ecuación.
El problema es que en este caso, no se puede añadir la variable
SlideName como factor fijo en la matriz de diseño ya que el
factor SlideNameSaraWTATMA6 no es estimable. Este aviso nos
indica que esta columna es redundante y puede expresarse como la
combinación de otras.
Por este motivo se emplea SlideName como bloqueo, lo que
permite controlar la variación técnica por slide.
## [1] "grupoIba1_Control_Primary" "grupoIba1_Control_Recurrence"
## [3] "grupoIba1_Nivolumab_Primary" "grupoIba1_Nivolumab_Recurrence"
## [5] "grupoSox2_Control_Primary" "grupoSox2_Control_Recurrence"
## [7] "grupoSox2_Nivolumab_Primary" "grupoSox2_Nivolumab_Recurrence"Se redefinen los nombres.
## [1] "Iba1_Control_Primary" "Iba1_Control_Recurrence"
## [3] "Iba1_Nivolumab_Primary" "Iba1_Nivolumab_Recurrence"
## [5] "Sox2_Control_Primary" "Sox2_Control_Recurrence"
## [7] "Sox2_Nivolumab_Primary" "Sox2_Nivolumab_Recurrence"Se obtienen 8 niveles de la variable grupo.
Al pasar el consensus.correlation a lmFit,
se “bloquea” el efecto técnico. Este evita que el modelo confunda la
variabilidad entre portaobjetos con diferencias biológicas reales,
aumentando mucho la potencia estadística para detectar cambios reales
entre Iba1 y Sox2.
corfit <- duplicateCorrelation(
v_limma,
design = design_limma,
block = colData(datos_lrb)$SlideName
)
fit_limma <- lmFit(
v_limma,
design_limma,
block = colData(datos_lrb)$SlideName,
correlation = corfit$consensus
)Se construye la matriz de contraste para el análisis de expresión diferencial (DE) entre Iba1 y Sox2 en situaciones de tumor primario y sin aplicar el tratamiento.
contr_matrix_limma <- makeContrasts(
Iba1_vs_Sox2_Control_Primary = Iba1_Control_Primary - Sox2_Control_Primary,
levels = colnames(design_limma)
)
fit_contrast_limma <- contrasts.fit(fit_limma, contrasts = contr_matrix_limma)
efit_limma <- eBayes(fit_contrast_limma, robust = TRUE)Se genera un volcano plot, para la mejor observación de estos resultados.
resultado_limma <- topTable(efit_limma, coef = "Iba1_vs_Sox2_Control_Primary", number = Inf, sort.by = "P")
resultado_limma$significant <- with(resultado_limma, adj.P.Val < 0.05 & abs(logFC) > 1)
ggplot(resultado_limma, aes(x = logFC, y = -log10(adj.P.Val))) + geom_point(aes(color = significant), alpha = 0.6) +
scale_color_manual(values = c("grey", "red")) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
theme_minimal() +
labs(
title = "Iba1+ vs. Sox2+ AOIs en pacientes Control con tumores primarios",
x = "log2 Fold Change",
y = "-log10(FDR)"
)
El Volcano Plot resultante del contraste Iba1+ vs. Sox2+ muestra una
asimetría hacia el enriquecimiento de genes en la población mieloide
(Iba1), con una significancia estadística extremadamente elevada (FDR
< \(10^{-40}\)). Esta distribución
confirma que la estrategia de bloqueo por SlideName ha
preservado las identidades biológicas de los segmentos, permitiendo que
emerjan con claridad tanto las firmas de presentación antigénica en
Iba1+ como los factores de progresión tumoral (p. ej., VEGFA) en
Sox2+.
7. GSEA y visualización con vissE
El Análisis de Enriquecimiento de Conjuntos de Genes (GSEA) es un método bioinformático que permite determinar si un conjunto de genes de un dataset está relacionado y muestra cambios de forma coordinada en una condición concreta.
En este pipeline se hará uso el método de GSEA del paquete limma que
utiliza la función fry().
Es necesario seleccionar colecciones de genes, donde se hayan definido previamente asociaciones entre ellos según su función biológica, su implicación en rutas metabólicas o vías de señalización.
Se seleccionan los siguientes conjuntos de genes para realizar un análisis de enriquecimiento de conjuntos de genes:
MSigDB Hallmarks: conjuntos de genes de la colección hallmarks de MSigDB.
MSigDB C2: conjuntos de genes de la colección C2 de MSigDB, que contiene conjuntos de genes curados, como aquellos obtenidos de bases de datos como BioCarta, KEGG, PID y Reactome, así como de experimentos de perturbación química o genética.
GO BP: procesos biológicos de la base de datos de ontología génica (gene ontology).
GO MF: funciones moleculares de la base de datos de ontología génica.
GO CC: componentes celulares de la base de datos de ontología génica.
FDR < 0.05 indica que el conjunto de genes está significativamente enriquecido.
7.1 Cargar conjunto de genes
Se cargan los conjuntos de genes utilizando el paquete
msigdb y se extraen únicamente los conjuntos de genes que
se han descrito.
library(msigdb)
library(GSEABase)
msigdb_hs <- getMsigdb(version = '7.2')
msigdb_hs <- appendKEGG(msigdb_hs) #Añadimos información de vías KEGG
sc <- listSubCollections(msigdb_hs) #Lista todas las subcolecciones dentro de MSigDB
gsc <- c(subsetCollection(msigdb_hs, c('h')),
subsetCollection(msigdb_hs, 'c2', sc[grepl("^CP:",sc)]),
subsetCollection(msigdb_hs, 'c5', sc[grepl("^GO:",sc)])) %>%
GeneSetCollection() #filtra las subcolecciones de interés y las organiza en un objeto con el formato estándar para análisis de enriquecimiento con funciones como fry o GSEA.7.2 Analisis de enriquecimiento
Se preprocesan los datos de los conjuntos de genes cargados,
filtrando los que contienen menos de 5 genes y creando una lista de
índices para formatear los resultados antes de aplicar
fry().
Se elige trabajar con el objeto v_limma correspondiente
al contraste Iba1 vs Sox2 en tumores primarios de pacientes control, con
el fin de caracterizar diferencias funcionales basales entre
compartimentos celulares.
fry_indices_limma <- ids2indices(lapply(gsc, geneIds), rownames(v_limma), remove.empty = FALSE)
names(fry_indices_limma) <- sapply(gsc, setName)
gsc_category_limma <- sapply(gsc, function(x) bcCategory(collectionType(x)))
gsc_category_limma <- gsc_category_limma[sapply(fry_indices_limma, length) > 5]
gsc_subcategory_limma <- sapply(gsc, function(x) bcSubCategory(collectionType(x)))
gsc_subcategory_limma <- gsc_subcategory_limma[sapply(fry_indices_limma, length) > 5]
fry_indices_limma <- fry_indices_limma[sapply(fry_indices_limma, length) > 5]
names(gsc_category_limma) = names(gsc_subcategory_limma) = names(fry_indices_limma)Ahora se aplica fry() con todos los conjuntos de genes
filtrados.
Se dividen los conjuntos de genes por categoría (gsc_category) y se aplica la función para cada categoría.
Se diseña una función para el output donde se combinan los resultados en un dataframe.
fry_indices_cat_limma <- split(fry_indices_limma, gsc_category_limma[names(fry_indices_limma)])
fry_res_out_limma <- lapply(fry_indices_cat_limma, function (x) {
limma::fry(v_limma, index = x, design = design_limma, contrast = contr_matrix_limma[,1], robust = TRUE)
})
post_fry_format_limma <- function(fry_output_limma, gsc_category_limma, gsc_subcategory_limma){
names(fry_output_limma) <- NULL
fry_output_limma <- do.call(rbind, fry_output_limma)
fry_output_limma$GenesetName <- rownames(fry_output_limma)
fry_output_limma$GenesetCat <- gsc_category_limma[rownames(fry_output_limma)]
fry_output_limma$GenesetSubCat <- gsc_subcategory_limma[rownames(fry_output_limma)]
return(fry_output_limma)
}
fry_res_sig_limma <- post_fry_format_limma(fry_res_out_limma, gsc_category_limma, gsc_subcategory_limma) %>%
as.data.frame() %>%
filter(FDR < 0.05)
head(fry_res_sig_limma)Se grafican los 10 conjuntos de genes Up y Down - Regulated en Iba1 frente Sox2. El conjunto de genes que se encuentran Up-regulated están más activos en Iba1, en cambio, los que se encuentran Down-regulated se están más activos en las zonas Sox2.
fry_res_sig_limma %>%
arrange(FDR) %>%
filter(Direction == "Up") %>%
.[seq(10),] %>%
mutate(GenesetName = substr(GenesetName, 1, 50),
GenesetName = factor(GenesetName, levels = GenesetName[order(-log10(FDR))])) %>%
ggplot(aes(GenesetName, -log(FDR))) +
geom_bar(stat = "identity", fill = "red") +
theme_bw() +
coord_flip() +
ggtitle("Up-regulated Iba1+ Condition")
fry_res_sig_limma %>%
arrange(FDR) %>%
filter(Direction == "Down") %>%
.[seq(10),] %>%
mutate(GenesetName = substr(GenesetName, 1, 50),
GenesetName = factor(GenesetName, levels = GenesetName[order(-log10(FDR))])) %>%
ggplot(aes(GenesetName, -log(FDR))) +
geom_bar(stat = "identity", fill = "blue") +
theme_bw() +
coord_flip() +
ggtitle("Down-regulated Iba1+ Condition")
Los conjuntos sobreexpresadsos en Iba1 muestran una significación estadística extrema (\(-log(FDR) > 100\)), lo que indica una señal biológica masiva y coherente con el linaje de macrófagos/microglía. Se destaca la fagocitosis de los patógenos, que es el conjunto de genes más robusto y que en un tumor indica que la microglía está intentando limpiar el entorno tumoral. La red Tyrobp que describe una activación profunda de las células y la regulación del superóxido (ROS), que implica actividad efectora pro-inflamatoria.
Los conjuntos de genes down-regulated son los que se sobreexpresan en Sox2 indicadores de células tumorales. Es coherente lo que se observa, ya que se asocian a una anomalía clínica donde se modifican las estructuras del tejido cerebral, también a la proliferación celular maligna y a genes en estados menos diferenciados.
Es un análisis correcto porque coincide con los resultados observados en el volcano plot. Indagando en las funciones de las vías, los genes “Up” definen funciones efectoras inmunes (Iba1), mientras que los genes “Down” definen la identidad del linaje neural y la progresión tumoral (Sox2). Esto confirma que se logra capturar con éxito la interacción entre el tumor y la respuesta inmunitaria del paciente.
En conclusión, el uso de GeoMx DSP combinado con un riguroso análisis bioinformático permite diseccionar la realidad de estas poblaciones celulares que coexisten en estrecha proximidad. Este enfoque hace posible perfilar sus funciones específicas, capturando con éxito la compleja interacción entre la progresión del linaje tumoral (Sox2+) y la respuesta inmunitaria del paciente (Iba1+).
8. Deconvolución celular
La Deconvolución Celular es una técnica que permite estimar las proporciones de los diferentes tipos celulares que se encuentran en una muestra de tejido.
En el paquete standR, se usa la función
prepareSpatialDecon para comunicar un objeto
SpatialExperiment al paquete SpatialDecon, con
el fin de realizar la deconvolución celular.
Sin embargo, SpatialDecon requiere sondas negativas para
establecer el ruido de fondo de los datos. El ruido de fondo es la señal
medida por el equipo que no proviene de la expresión real de los genes,
sino de artefactos técnicos.
Por lo tanto, es necesario reconstruir el objeto
SpatialExperiment usando el parámetro
rmNegProbe = FALSE, lo que desactiva la eliminación de las
sondas negativas, y luego volver a ejecutar los pasos de control de
calidad (QC).
Se realiza la preparación de los datos, desde un punto de partida muy
inicial: leer datos, QC por ROI, QC conteo núcleos >150,
normalización por TMM. En este caso se puede utilizar directamente los
datos_tmm.
library(SpatialDecon)
colData(datos_tmm)$grupo <- paste(colData(datos_tmm)$SegmentLabel, colData(datos_tmm)$trial_setting, colData(datos_tmm)$tumor_setting, sep = "_")
datos_deco <- prepareSpatialDecon(datos_tmm)El objeto que se obtiene con la función
prepareSpatialDecon() contiene dos matrices, la matriz de
conteos normalizados y la matriz del ruido de fondo.
## [1] "normCount" "backGround"A continuación, se prosigue con la ejecución de
SpatialDecon, el paquete de R cuyo algoritmo matemático
permite estimar las proporciones celulares gestionando simultáneamente
el ruido de fondo.
Para este análisis se utiliza la matriz llamada SafeTME
como base de referencia, que está diseñada para estimar células inmunes
y del estroma en el microambiente tumoral. Las matrices son diccionarios
que relacionan el tipo y la cantidad de cada gen que están expresados en
cada linaje celular. La matriz Safe_TME ha sido creada
específicamente para evitar genes que se expresan comúnmente en células
cancerosas, de modo que la estimación de los tipos celulares no se vea
sesgada.
Mediante un modelo de regresión lineal, el algoritmo utiliza esta referencia para desglosar la señal total de cada AOI, calculando qué proporción de cada tipo celular explica de forma óptima el perfil de expresión observado en la muestra.
data("safeTME")
heatmap(sweep(safeTME, 1, apply(safeTME, 1, max), "/"),
labRow = NA, margins = c(10, 5))
Se aplica la función spatialdecon() para los datos
normalizados, teniendo en cuenta el ruido de fondo y el catálogo
SafeTME.
res <- spatialdecon(norm = datos_deco$normCount,
bg = datos_deco$backGround,
X = safeTME,
align_genes = TRUE)Se crean subconjuntos de las zonas celulares Iba1 y Sox2, para una observación más ordenada que permite una mejor comparación.
samples_subset_iba1 <- colnames(datos_tmm)[colData(datos_tmm)$SegmentLabel == "Iba1"]
subset_prop_iba1 <- res$prop_of_all[,samples_subset_iba1]
spe_sub_iba1 <- datos_tmm[,samples_subset_iba1]
samples_subset_sox2 <- colnames(datos_tmm)[colData(datos_tmm)$SegmentLabel == "Sox2"]
subset_prop_sox2 <- res$prop_of_all[,samples_subset_sox2]
spe_sub_sox2 <- datos_tmm[,samples_subset_sox2]Se crean los gráficos.
library(tidyverse)
library(dplyr)
# 1. Convertir la matriz de proporciones a formato largo
df_plot_iba1 <- as.data.frame(t(subset_prop_iba1)) %>%
rownames_to_column("AOI") %>%
pivot_longer(cols = -AOI, names_to = "CellType", values_to = "Proportion")
# 2. Extraer la información de los grupos desde el colData
grupos_iba1 <- as.data.frame(colData(spe_sub_iba1)) %>%
rownames_to_column("AOI") %>%
dplyr::select(AOI, grupo) # Ajusta 'Grupo' al nombre real de tu columna
# 3. Unir las proporciones con los grupos
df_final_iba1 <- left_join(df_plot_iba1, grupos_iba1, by = "AOI")
ggplot(df_final_iba1, aes(x = AOI, y = Proportion, fill = CellType)) +
geom_bar(stat = "identity", width = 1) +
facet_wrap(~grupo, scales = "free_x", ncol = 4) +
theme_minimal() +
labs(x = "Muestras por Grupo",
y = "Proporción Celular",
fill = "Tipo Celular") +
theme(
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.spacing = unit(0.1, "lines"), # Une un poco más las barras
strip.background = element_rect(fill = "gray90"), # Estilo para la etiqueta del grupo
legend.position = "bottom"
)
# 1. Convertir la matriz de proporciones a formato largo
df_plot_sox2 <- as.data.frame(t(subset_prop_sox2)) %>%
rownames_to_column("AOI") %>%
pivot_longer(cols = -AOI, names_to = "CellType", values_to = "Proportion")
# 2. Extraer la información de los grupos desde el colData
grupos_sox2 <- as.data.frame(colData(spe_sub_sox2)) %>%
rownames_to_column("AOI") %>%
dplyr::select(AOI, grupo) # Ajusta 'Grupo' al nombre real de tu columna
# 3. Unir las proporciones con los grupos
df_final_sox2 <- left_join(df_plot_sox2, grupos_sox2, by = "AOI")
ggplot(df_final_sox2, aes(x = AOI, y = Proportion, fill = CellType)) +
geom_bar(stat = "identity", width = 1) +
facet_wrap(~grupo, scales = "free_x", ncol = 4) +
theme_minimal() +
labs(x = "Muestras por Grupo",
y = "Proporción Celular",
fill = "Tipo Celular") +
theme(
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.spacing = unit(0.1, "lines"), # Une un poco más las barras
strip.background = element_rect(fill = "gray90"), # Estilo para la etiqueta del grupo
legend.position = "bottom"
)
Los graficos se dividen en paneles según los 8 niveles posibles de la
variable grupo. Cada barra vertical corresponde a una AOI y
muestra la proporción de cada tipo celular mediante una escala de
colores.
Se observa una distinción clara de las proporciones celulares entre las dos zonas estudiadas. En las AOI de Iba1 se encuentran mayoritariamente proporciones de macrófagos, confirmando así la identidad mieloide de estos segmentos. Mientras que en las AOI de Sox2 hay una proporción muy alta de células T (CD4+ y CD8+), este hallazgo representa que las células cancerígenas han sido infiltradas por los linfocitos. a pesar de eso, no se observan diferencias significativas entre pacientes tratados y o tratados, ni en la deconvolución ni como se ha mostrado anteriormente en el análisis diferencial, lo que se podría clasificar como una presunta infiltración ineficaz.
Se observa una mayor variabilidad intragrupo en pacientes recurrentes que en pacientes primarios, lo que sugiere una mayor heterogeneidad celular tras la recurrencia del tumor.
Se observa, también, una capa de células endoteliales y fibroblastos constante en ambas zonas.
9. Análisis de proporción diferencial
El análisis de proporción diferencial tiene como objetivo identificar
estadísticamente cambios en la abundancia relativa de tipos celulares
entre diferentes condiciones experimentales. Para ello, se emplea la
librería speckle, que implementa la función
propeller().
En este apartado se analiza la distinción entre los nichos Iba1 y
Sox2 en pacientes control con tumores primarios. Los resultados de la
deconvolución res se procesan mediante la función
convertDataToList() para adaptar las proporciones a un
formato que sea compatible con propeller.ttest().
El fundamento matemático de este test reside en la transformación de
las proporciones mediante la función arcoseno (asin), la
cual estabiliza la varianza y permite que el test sea más robusto.
Finalmente, se aplica el contraste entre las poblaciones celulares para
identificar aquellas poblaciones celulares que varían significativamente
(FDR < 0.05) entre ambos nichos.
library(speckle)
datos_tmm_df <- datos_tmm
samples_full <- colnames(datos_tmm_df)
prop_full <- res$prop_of_all[, samples_full]
spe_full <- datos_tmm[, samples_full]
propslist <- convertDataToList(prop_full,
data.type = c("proportions"),
transform="asin",
scale.fac=colData(spe_full)$AOINucleiCount)
design_prop <- model.matrix(~ 0 + grupo, data = as.data.frame(colData(spe_full)))
colnames(design_prop) <- gsub("^grupo","",colnames(design_prop))
contr <- makeContrasts(
Iba1_vs_Sox2_Control_Primary = Iba1_Control_Primary - Sox2_Control_Primary,
levels = colnames(design_prop)
)
outs <- propeller.ttest(propslist, design_prop, contr, robust=TRUE,trend=FALSE, sort=TRUE)Se realiza la visualización de las proporciones mediante violin plots.
diff_ct <- outs %>%
filter(FDR < 0.05) %>%
rownames()
colData(spe_full)$samples_id <- rownames(colData(spe_full))
prop_full[diff_ct,] %>%
as.data.frame() %>%
rownames_to_column("CellTypes") %>%
gather(samples, prop, -CellTypes) %>%
left_join(as.data.frame(colData(spe_full)), by = c("samples"="samples_id")) %>%
ggplot(aes(SegmentLabel, prop, fill = CellTypes)) +
geom_violin() +
facet_wrap(~CellTypes) +
theme_bw() +
xlab("Iba1 vs Sox2 en pacientes control con tumores primarios") +
ylab("Proportion")
Los resultados muestran las diferencias en la composición celular entre las zonas Iba1 y Sox2 para el grupo de pacientes control con tumores primarios. Tras el análisis, destacan dos hallazgos principales con significación estadística:
Macrófagos: Presentan un enriquecimiento marcado en las zonas Iba1, lo que confirma la eficacia del marcador para aislar nichos mieloides.
Células T CD8+: Se encuentran en una proporción significativamente mayor en las zonas Sox2, corroborando una vez más la infiltración de linfocitos citotóxicos dentro del compartimento tumoral.
Por el contrario, el resto de los linajes celulares (como fibroblastos o células endoteliales) mantienen proporciones similares en ambas regiones.
10. Conclusiones
El análisis integral de los datos de transcriptómica espacial de GeoMx DSP mediante el flujo de trabajo de GeoMxTools y StandR, permite extraer las siguientes conclusiones fundamentales:
- Robustez del Pipeline Bioinformático:
El pipeline cuenta con métodos robustos de preprocesamiento, como la estrategia de normalización y corrección de efectos técnicos, que permiten reducir la variabilidad técnica y potenciar las diferencias biológicas entre condiciones.
- Validación de la Identidad de los segmentos celulares (Iba1 y Sox2):
Se ha confirmado mediante los análisis posteriores (expresión diferencial, análisis de enriquecimiento, deconvolución celular y análisis de proporción diferencial) que la selección de segmentos fue precisa.
- Evidencia de Infiltración Ineficaz:
Uno de los hallazgos más relevantes del análisis de este dataset es la detección de una elevada proporción de linfocitos T (CD4+ y CD8+) dentro de las zonas tumorales (Sox2+). A pesar de esta infiltración, la ausencia de diferencias transcriptómicas significativas entre pacientes tratados con Nivolumab y controles sugiere un estado de infiltración funcionalmente ineficaz. Este hallazgo coincide con el estudio publicado de este dataset.